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Abstract. The dynamical Casimir effect for a massless scalar field in 1+f -dimensions 
is studied numerically by solving a system of coupled first-order differential equations. 
The number of scalar particles created from vacuum is given by the solutions to this 
system which can be found by means of standard numerics. The formalism already used 
in a former work is derived in detail and is applied to resonant as well as off-resonant 
cavity oscillations. 



PACS numbers: ll.lO.-z, 42.50.Pq 
1. Introduction 

The possibility of creating photons out of vacuum fluctuations of the quantized 
electromagnetic field in dynamical cavities, the so-called dynamical Casimir effect (see 
^ for a review), demonstrates the highly non-trivial nature of the quantum vacuum. 
A scenario of particular interest are so-called vibrating cavities 2^ where the 



distance between two parallel (ideal) mirrors changes periodically in time. The 
occurrence of resonance effects between the mechanical motion of the mirror and the 
quantum vacuum leading to (even exponentially) increasing occupation numbers in 
the resonance modes makes this configuration the most promising candidate for an 
experimental verification of this pure quantum effect. 

Particle creation in one- dimensional vibrating cavities has been studied in numerous 
works 13 in El El 13 IHl El mH HI]- When considering small amplitude oscillations, 
analytical results can be deduced showing that under resonance conditions the particle 
occupation numbers increase quadratically in time. Particle creation due to off-resonant 
wall motions has been investigated in, e.g., 0. The evolution of the energy density in a 
one-dimensional cavity with one vibrating wall has also been studied by many authors 
[121 Uni Hn Uni Uni mi UHI demonstrating that the total energy inside a resonantly 
vibrating cavity grows exponentially in time (see also |^) while the total particle 
number increases only quadratically. Thus a pumping of energy into higher frequency 
modes takes place and particles of frequencies exceeding the mechanical frequency of 
the oscillating mirror are created. The energy for this process is provided by the energy 
which has to be given to the system from outside to maintain the motion of the mirror 
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against the radiation reaction force 1201 122] • The more reahstic case of a three- 
dimensional cavity is studied in|221l2Sll21l23l2nil2Ill2Hll2ni- Field quantization inside 
cavities with non-perfect boundary conditions has been investigated in, e.g., [SOI EI] and 
corrections due to finite temperature effects are treated in [221 EBl IM] • The question of 
how the quantum vacuum interacts with the (classical) dynamics of the cavity has been 
addressed in |I2l EE EH ■ 

In this work we present a formalism allowing for numerical investigation of the 
dynamical Casimir effect for scalar particles in a one-dimensional cavity. (For related 
numerical work see also pHj EHl EH]-) We introduce a particular parametrization for 
the time evolution of the field modes yielding a system of coupled first-order differential 
equations. The solutions to this system determine the number of created particles and 
can be obtained by means of standard numerics. We employ the formalism to investigate 
the creation of real massless scalar particles in a resonantly as well as off-resonantly 
vibrating cavity and compare the numerical results with analytical predictions. These 
results are complementary to the ones already presented in |3T]. 

With this formalism at hand the dynamical Casimir effect can be investigated fully 
numerically making it possible to study a variety of scenarios where no analytical results 
are known (large amplitude oscillations, arbitrary wall motions etc.). Of special interest 
is of course the realistic case of the electromagnetic field in a three-dimensional cavity. 
Being easily extendable to arbitrary space dimensions the presented formalism can be 
used also in this case. In particular it allows to calculate numerically the TE-mode 
contribution [2ZI to the photon creation taking the influence of the intermode coupling 
fully into account jl21- Hence the formalism can be used to cross-check analytical results 
also in this realistic case which might be of importance for future experiments. Let us 
finally note that the tools (and the numerical formalism presented here) used to study 
the dynamical Casimir effect can also be employed to investigate graviton generation in 
braneworld cosmology 



2. Hamiltonian and equations of motion 

We consider the Hamilton operator 

= \Y.[pI + - ^ E Mnrr.{t) [qnPm + PrnQn] (l) 

with 

Pn = Qn + y^QniMmnit) and Mnm{t) = / c/x 0„(t, x)0^(t, x) (2) 

Jl(t) 

describing the dynamics of a massless real scalar field $ = $(t, x) on a time-dependent 
interval I{t) = [0, l{t)] in terms of the canonical operators and Pn subject to the usual 
equal-time commutation relations. The corresponding canonical variables g„ and p„ are 
introduced via the expansion $(t,x) = J2n (lnit)(pn{t,x) of the field and its momentum 



Numerical approach to the dynamical Casimir effect 



3 



n(t,x) = Y,n Vnif)4'nifix) in time-dependent (instantaneous) eigenf unctions (j)n(t,x) 
satisfying the eigenvalue equation 

-dlMf^) = ^lit)Mt,x) (3) 

on I(t) with time- dependent eigenvalues fin{t) [7j. The overdot denotes the derivative 
with respect to time t and we are using units with h = c = 1. The time-dependent 
so-called coupling matrix Mnm arises due to the time- dependent boundary condition for 
the field at x = l{t) enforcing the eigenfunctions of —d^ to be time-dependent. We take 
the boundary conditions at x = and x = l{t) to be of the form 

[ai* + a2d,<^] U=o = + b2d,^] \,=i^t) = V t (4) 

with constants ai,a2,&i and 62 ensuring that the set of eigenfunctions {(j)nit,x)} is 
complete and orthonormal for all times. 

Adopting the Heisenberg picture, the equations of motion for Qnit) read lUIZI 

g„ + Qlqn + [Mmn - Mnm] Qm + Y^ [^mn " Nnm = (5) 

m m 

where Nnm = J2k MnkMmk- The structure of the intermode coupling mediated by the 
coupling matrix Mnm depends on the particular kind of boundary conditions which 
decide on the specific form of the instantaneous eigenfunctions x). It is worth 
noting that the Hamiltonian does not correspond to the energy of the field because 
the coupling term does not contribute to the total energy defined via the energy 
momentum tensor |7 . From the Hamiltonian and the equations of motion (0) 
we identify two external time dependences in the equations which will lead to particle 
creation: (i) the time-dependent eigenfrequencies and (ii) the coupling matrix 

Mnmit), called the squeezing and acceleration effect, respectively. 

3. Vacuum and particle definition 

Let us assume that the motion of the wall is switched on at t = with l{t) following 
a prescribed trajectory for a duration ti, ceases afterwards and is at rest again. Before 
and after the motion the coupling matrix vanishes and the time evolution of the 
operator g„ is determined by the equation of an harmonic oscillator with constant 
frequency = f2„(t < 0) and = > ti), respectively |. The corresponding 

Hamilton operator describing the quantized field for t < and t > ti can then 
be diagonalized by introducing time-independent annihilation and creation operators 
{an, all}, corresponding to the particle notion for t < 0, and {y4„,y4|^} associated with 
the particle notion for t >ti via § 

t<0: qnit) = + h.c, Pn{t) = t\-^ai e^^°* + h.c, (6) 

I Here the final position l(ti) = li is assumed to be arbitrary. In case of a vibrating cavity it is natural 
to consider times ti after which the dynamical wall has returned to its initial position. 
§ We are assuming that il^^ ^ and il^ ^ for all n. 
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t > : Ut) = " , + h.c, pnit) = ^\H^Al e^^-^*'*^) + h.c..(7) 



The initial and final vacuum states |0,t < 0) and \0,t > ti), respectively, are introduced 
as the ground states of the corresponding diagonal Hamilton operators: 

f [alon + 1/2] with |0, t < 0) = for t < 

H = E\ ■ (8) 

" [ nj^ [AlAn + 1/2] with An |0, t > ti) = for t > ti 

The set of initial-state operators {a„, dl} is related to the set of final-state operators 
{An, Al^} by a Bogoliubov transformation 

An = Y^ [Amnih) d^ + I3^n{tl) (9) 

m 

where Amn(ti) and Bmniti) satisfy the relations 

['^mn-^mA: ~ ^mn^mk] = "^nfc ; ^ {A-mn^mk ~ ^mn-^mk] = 0. (10) 

m m 

For t > ti the particle number operator iV„ = Aj^An defined with respect to the final 
vacuum state counts the number of physical particles. The number of particles created 
in a mode n during the motion of the wall is given as the expectation value of Nn with 
respect to the initial vacuum state \0,t < 0): 

Nnih) = {0,t< 0\Nn\0,t< 0)=Y^ \Bmnitl)\^ (11) 

m 

Accordingly, total particle number N{ti) and energy E{ti) of the motion induced 
radiation are given by 

iV(ti) = ' ^(^l) = E^n^"(^l)- (12) 

n n 

Both quantities are in general ill defined and therefore require appropriate regularization. 
For a time dependence of the boundary l{t) which is not sufficiently smooth, i.e. it 
exhibits discontinuities in its time- derivative appearing for instance when switching the 
motion on and off instantaneously, one may expect that part of the particle creation 
is due to this discontinuity in the velocity which may cause the excitation of modes of 
even arbitrary high frequencies. Hence a (large) contribution to the predicted particle 
creation may be spurious and in the case that arbitrary high frequency modes become 
excited the summations in (fT^ do not converge. This can be avoided most easily by 
introducing a frequency cut-off which effectively smoothes the dynamics l{t). When 
calculating the quantities (|T2|l numerically we will make use of such a frequency cut- 
off which is determined by the stability of the numerical results for single modes, i.e. 
stability of the expectation value ()11|) with respect to the cut-off. Note that an explicit 
frequency cut-off also accounts for imperfect (non-ideal) boundary conditions for high 
frequency modes 
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4. Time evolution 



During the motion of tlie boundary some or even infinitely many modes may be coupled. 
For t > the operators qn{t) and p„(t) given by U^qn{0)tj and U'^Pn{0)U, respectively, 
with U = ij{t,0) = Texp(— z /q dt'H{t')) and T denoting the time-ordering operator, 
can be expanded in initial state operators a„, a]^ and complex functions e^^\t): 



Ut > 0) = uHnio)u = ^ 

m 

Pn{t > 0) = U^Pn{0)U = J2 



ei™)(t)+h.c.. 



2Dt 



(m) , 



(13) 
+ h.c(14) 



By using the Heisenberg equation 0(t) = i[H{t),0{t)] + {dO{t) / dt)e^^\, it is 
straightforward to show that the functions e^™) {t) satisfy the same differential equation 
(jni) as qnif)- Notice that insertion of Eq. (fT!?|l into the mode expansion for $ leads to 
the decomposition of the field used in, e.g., [H El IH|- Through the formal expansion 
(fT!?|l we have reduced the problem of finding the time evolution for the operator qnif) 
to the problem of solving the system of coupled second-order differential equations (0) 
for e^^\t). Demanding that Eqs. (fT^ and p4|) have to match with the corresponding 
expressions (jHl) at t = leads to the initial conditions 



^(m) 



(0) = (5„ 



(0) 



(15) 



Hence with Mm„(0) vanishing only if /(O) = the initial condition e^'")(0) is not simply 



—iVt^^5nm when dealing with boundary motions l{t) which have a discontinuity in the 
velocity at t = 0. Matching ((Zj) with (fT!?|) and (fT^ at t = ti one finds 

1 



n 



_(m) 



no p" ni 



■Am) 



(m), 



(16) 



(17) 



Starting with the initial vacuum |0,t < 0) the Bogoliubov transformation © has to 
become trivial for ti = 0, i.e. = a^, implying the vacuum initial conditions 

Amnio) = 6„,n and Bmn{0)=0 (18) 

which are consistent with the initial conditions (jl5|) . The emergence of Mmn{0) in the 
initial conditions (fT^ therefore guarantees to meet the vacuum initial conditions when 
the motion of the boundary starts instantaneously with a non-zero velocity. 
By introducing the auxiliary functions || 



,(m) 



,(m) 



n 



r>(«)+5:A4„(«)ei:"'(«) 



{t) + Y,MUt) 

k 



(ra) 



it) 



(19) 
(20) 



A derivation can be found in Appendix A. 
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the expressions and (|17|) can be rewritten as 

~^^:-\t,)+A-{t^)vt\t, 



with 



•^mn (tl 



1 ± 



(21) 

(22) 
(23) 



The quantity A^(ti) is a measure for the deviation of the final state of the cavity, 
characterized by the cavity length /(ti), with respect to its initial size Iq. If at time 
tl the cavity size is equal to the initial size lo, for instance in the important case 
that tl is a multiple of the period of oscillations of the cavity, we have Bmn{ti) = 
{l/2)^nynU:-\ti) and therefore = {l/^)EM^lWrr\ti)\'. 

The advantage of introducing the functions and r]^"^-* is that they satisfy the 
following system of first-order differential equations: 



m) 



(m) 



with 



1 ± 



(m) 



nut) 



and c^fc(t) 



nk 



it) 



(24) 
(25) 

(26) 



Besides the time-dependent frequency Qn(t) only the coupling matrix Mknif) enters into 
this system of coupled differential equations but neither Nnkif) nor its time derivative 
Mkn{t). The vacuum initial conditions (fTH|) entail the initial conditions for the functions 
^^"'^ and to be 

e^Ho) 



25 r. 



0. 



(27) 



Let us stress that all derivations and equations shown so far, do not rely on particular 
symmetry properties of the coupling matrix. 

By means of Eq. (j22p the number of particles created from vacuum during the 
dynamics of the cavity as well as the associated energy may now be calculated from the 
solutions ^^'"^ and t]^"*-* of the system of coupled first-order differential equations formed 
by Eqs. (j^ and (j^ . 

In order to obtain the numerical results presented in the next section we proceed 
in the following way: A cut-off parameter fcmax is introduced to make the system of 
differential equations finite and suitable for a numerical treatment. The system of 
coupled differential equations is then evolved numerically from t = up to a final time 
tmax and the expectation value (fTT|) is calculated for several times in between. By doing 
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so we interpret ti as a continuous variable such that Eq. (fTTj) becomes a continuous 
function of time^. Consequently, the stability of the numerical solutions with respect to 
the cut-off has to be ensured. In particular /cmax will be chosen such that the numerical 
results for the number of particles created in single modes are stable. Furthermore, 
the quality of the numerical results is assessed by testing the Bogoliubov relations (fTnjl . 

This procedure is of course not without problems when the expectation values are 
evaluated also for times ti at which l{ti) 7^ 0. The used particle definition requires then 
a matching of the solutions to expressions corresponding to the static configuration with 
l(ti) = 0, hence a discontinuity in the velocity appears which may cause spurious effects. 
However the cut-off automatically ensures that possible spurious effects do not yield a 
divergent total particle number (see also section 3). Indeed we will see that in the 
particular scenario of interest - vibrating cavity - the influence of this matching problem 
is tiny and the numerical results agree perfectly with analytical predictions. 

5. Numerical results and discussion 

In this section we consider a massless real scalar field subject to Dirichlet boundary 
conditions at x = 0, l(t) and the much studied sinusoidal cavity motion 



foin^m and Mnn{t) = with n,m = 1, 2, 3, ... . The motion whose absolute value 
of the velocity is maximal at the beginning of the motion as well as for times at which the 
wall returns to its initial position features the above described matching problem. In [H] 
we have already studied particle creation caused by this motion for the main resonance 
case u = with the same formalism. We have found that for sufficiently small e 
and appropriate fcmax the numerical results are in excellent agreement with analytical 
predictions of j3] and [3]. Furthermore, the influence of the initial discontinuity in the 
velocity of the motion (|28p has been investigated showing that it is negligible for e -C 1. 

Here we want to concentrate on higher resonances u = 2f2° with n > 1 and off- 
resonant frequencies (detuning). In the simulations we set /o = 1 and e = 0.001. For 
these parameters it is shown in [H] that the numerical results agree very well with 
analytical predictions derived under the assumption e ^ 1. The numerical results are 
compared with analytical expressions obtained in [HI IE] • Remarks about the numerics 
can be found in Appendix C. 

In Fig.H(a) we show the numerical results for the total particle number in the time 
range [0, 250] for resonant cavity frequencies u = 2n° = 2n'ir for n = 1.5, 2, 2.5 and 3 
and the associated energy of the created quantum radiation is depicted in Fign(b). The 

% Interpreting ti as a continuous function of time one can of course derive a corresponding system of 
coupled differential equations for Amn and Bmn (see Appendix B). 



l{t) = /o [1 + e sin(tut)] , e < 1. 
The time-dependent frequency and coupling matrix are given by [7] 



(28) 




(29) 
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(b) 




(a) 



— analytical 



n=3 
n=2.5 



n=1.5 



300 



Figure 1. (a) Total number of particles produced in a cavity vibrating with H28|l 
and uj ~ 2m: with n = 1.5,2,2.5 and 3. The small plot shows the results in 
the time range [0, 50] together with the analytical prediction (solid line) N{t) = 
n{4:n'^ - l)(10-37rt)7l2 of [i El valid for short times (IQ-^Tri) < 1. (b) Numerical 
results for the created energy E{t) corresponding to (a) together with the analytical 
prediction E{t) = [Ari^ - l)7r sinh^(nl0"37rt)/12 of 8 (solid line). The results 
correspond to the largest cut-off parameters as given in Fig. |5J 



corresponding particle spectra are shown in Fig. El for different cut-off parameters fcmax 
to demonstrate numerical stability of the results. Here stability of the numerical results 
means that for the lowest modes k the value Nk{t) remains unchanged (within numerical 
precision) under variation of fcmax- The spectra confirm that no modes k = 2np with 
p = 1,2, 3, ... are coupled (and therefore excited) as predicted by the coupling condition 
u = ± derived and discussed in [231 • 

For short times ent = 10~^7rt -C 1 the numerical results are well described by 
the analytical predictions of [HI E] • The numerically calculated spectra for times t = 25 
shown in Fig.|21are well fitted by the analytical expression Nk{t) = {2n — k)k{10~^7rtY/4 
for k < 2n and Nkit) = otherwise 0, predicting a parabolic shape of the 
particle spectrum. More quantitatively, for n = 2, for instance, the predicted values 
Ni{t = 25) = N-iit = 25) ~ 4.63 x lO^^, iV2(t = 25) = 6.17 x 10-=^ agree 
well with the values Ni{t = 25) = 4.62 x 10^^ N2{t = 25) = 6.14 x 10"^ and 
N^{t = 25) = 4.59 X 10^^ obtained from the numerical simulations with A;max = 50. 
Consequently, the total number of created particles is perfectly described by the 
expression N{t) = n(4n^ — l)(10~^7rt)^/12 [3 El as it is demonstrated in the small 
plot in Fig. m (a). 

For the entire integration range [0, 250] we compare the numerical results for 
the total energy associated with the created quantum radiation with the analytical 
expression E{t) = [An^ — l)7r sinh^(ra 10^'^7rt)/12 [HI predicting that the energy increases 
exponentially with time [Fig. E(b)]. The numerical values and the analytical prediction 
agree very well for n = 1.5 and 2. In the case of n = 2.5 and 3 we observe slight deviations 
towards the end of the integration range. This is due to the numerical instabilities in 
the corresponding particle spectra [cf Figs. 121(c) and (d)]. The numerical values for A^^ 
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k 



6 10 14 18 22 26 30 34 38 42 46 50 

k 





Figure 2. Particle spectra for (a) ui = Stt, (b) uj = Air, (c) lu = 5n and (d) uj = 
corresponding to the results shown in Fig. ^ The small plots compare the numerical 
results for Nk{t = 25) with the analytical prediction Nk{t) = {2n - fc)fc(10"^7rt)^/4 of 
[HI plotted for continuous values of k (solid line). 



with k larger than some value {k > 10 for n = 3, for instance) do not remain unchanged 
when varying /Cmax- Even A^^ is small for the higher frequencies compared to the values 
of for the excited lowest modes their contribution to the total energy is significant 
because of their high frequency. Hence relatively small instabilities in A^^^ for larger k 
give rise to a non-stable (with respect to fcmax) result for the energy. In order to gain 
better agreement of the numerical results for the energy for n = 2.5 and 3 with the 
analytical prediction a further increase of /cmax is necessary. 

We now consider the case of detuning uj = 27r(n + 6n). In an off- resonant 
vibrating one-dimensional cavity the total energy associated with the created particles 
may increase exponentially E{t) = (tt / 12) {An"^ - 1) sinh^(nVl - i^nlQ-^ t)/(l - 7^) 
if 7 < 1, quadratically E{t) = {n/3){An^ - l){mTlO-'H/2f if 7 = 1 or oscillate 
E{t) = (7r/12)(4n2 - 1) sm^{ny/i^ - IvrlQ-^ t)/{-^^ - 1) if 7 > 1, depending on the 
strength of detuning 6n parametrized by 7 = 5n x 10^/n jH]. In Fig. ini(a) results for the 
total energy obtained in simulations with different off-resonant frequencies are shown 
covering all three different possibilities for 7 and compared to the analytical predictions. 
In all cases the numerical results are very well described by the analytical expressions. 
Figure El (b) depicts the periods of the energy (and particle number) oscillations 
as obtained from the simulations and compares them with the analytical prediction 
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Figure 3. (a) Energy associated with the particles created in an off-resonantly 
vibrating one-dimensional cavity. Numerical results are compared to the analytical 
predictions of |S] • The numerical results are always shown for cut-off parameters fcmax 
which ensure numerical stability, (b) Period of particle number oscillations caused by 
detuning. The numerically obtained period is compared with the analytical prediction 
to = 10"^/ {ny^'j^ — 1) of 8 . (c) Maximal amplitude N{to/2) of the corresponding 
particle number oscillations fitted to the power law N{to/2) cx 



to = 10'^/ {riy/'-f'^ — 1) showing that both are in good agreement. The numerical values 
for the maximal amplitudes A^(to/2) of the corresponding particle number oscillations 
are shown in Fig. 121 (c) and fitted to the power law N(tQ/2) oc with values of a 

as indicated in the figure. 

The numerical results presented in this section are entirely in very good agreement 
with the corresponding analytical predictions derived for small amplitude oscillations 
e -C 1 which demonstrates the reliability of the numerical simulations. However, a few 
critical comments are in order. In our considerations above the analytical expressions 
have been treated as continuous functions of time. But strictly speaking, they are vahd 
only for times at which the moving wall has returned to its initial position. Moreover, in 
the numerical simulations the expectation values have been calculated also for times at 
which the velocity of the moving wall is non-zero (matching problem). Consequently one 
may expect that part of the particle production is spurious and in particular if modes 
of arbitrarily high frequency are excited, for instance due to the initial discontinuity 
in the velocity of the motion (j^ . and fcmax co the particle number diverges and 
the numerical results do not agree with the analytical predictions. In the numerical 
simulations this is automatically avoided due to the cut-off /cmax- Nevertheless, working 
with a finite cut-off is well motivated because it simulates imperfect boundary conditions 
for high frequency modes and, as a matter of course, is a necessity for a numerical 
treatment. For the numerical results presented above these spurious effects are negligibly 
small and therefore the numerical results agree very well with the analytical predictions. 
This is due to the fact that we restrict ourselves to small amplitudes e <^ 1 and hence to 
small velocities. Therefore the effect of the discontinuity in the velocity of the boundary 
motion on the particle creation is expected to be small. This has been studied for 
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the initial discontinuity in the velocity in jH] in detail demonstrating that the matching 
problem becomes only important for larger e. In Appendix C we discuss the convergence 
of the numerical results in comprehension with the numerical accuracy of the simulations 
for the case u; = Svr more detailed. 

6. Conclusion 

A formalism allowing for numerical investigation of particle creation from vacuum 
in dynamical cavities, i.e. the dynamical Casimir effect, has been presented. By 
introducing a particular parametrization for the time-evolution of the field modes inside 
the dynamical cavity, a system of coupled first-order linear differential equations has 
been derived. Physical quantities like the number of particles created during the 
dynamics of the cavity and the associated energy are determined by the solutions to 
this system which can be found by applying standard numerics. 

In continuation of the work ^] we have studied the creation of massless scalar 
particles due to resonant as well as off-resonant sinusoidal oscillations of one of the 
cavity walls. The numerical results are entirely in agreement with the analytical 
predictions derived in [HI E] demonstrating that the numerical simulations are reliable 
and the method introduced is appropriate to study the dynamical Casimir effect fully 
numerically. 

Potential problems inherent in the method, in particular the matching problem due 
to discontinuities in the velocity of the boundary motion yielding spurious contributions 
to the total particle number, have been discussed. It has been shown that this effect is 
negligibly small for cavity vibrations with a sufficiently small amplitude. 

Being derived very generally, the method is applicable for different kinds of 
boundary conditions of the form provided that the spectrum {Vtn} contains no 
zero mode, i.e. fi„ > Vra and can easily be extended to massive scalar fields by 
substituting for the frequency f2„ the corresponding expression for a massive scalar 
field, i.e. Vtn = \J {nny/P + m? where m is the mass. Furthermore, the generalization 
to higher dimensional cavities is straightforward. This makes it possible to study the 
dynamical Casimir effect for a variety of possible interesting scenarios where less or 
even nothing is known analytically. As already mentioned in the introduction, TE- 
mode photon creation in a three-dimensional rectangular cavity [23 can be studied 
with the same method as well because it can be related to the production of massive 
scalar particles in a one-dimensional cavity jl2]. However, more complicated boundary 
conditions than appearing for example when studying TM-mode photons ^Tj\ cannot 
be treated within this approach. 
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Appendix A. Derivation of ^'f^'> and rj'^'^ 

The auxiliary functions ^i"'^^) and 77l™)(t) [Eqs. (HH)), ^] can be introduced in 
the following way: Define an operator 6„(t) via hn{t) := f/^(t, 0) a„ f/(t, 0) with a„ 
being the annihilation operator corresponding to the initial state [Eq. ©], i.e. a„ = 
(l/v^)[ynOg„(0) + (Vv^)Pn(O)]. Hence kit) = (l/v^)[yfiOf/t(t, 0)g„(0)f/(t, 0) + 

(i/w^)[/t(t,0)p„(0)f/(t,0)]. By using Eqs. (O and ^ one derives 



kit > 0) = 5] 



^2. 



(A.i; 



with ^^rTKt) and r]^^\t) defined in Eqs. (0 and Note that this definition of 

the time evolution for does not account for an explicit time-dependence of bnit). 
Therefore, in general, no meaningful notion of particles may be assigned to the operator 
bnit). This manifests itself in the relation between the operator An corresponding to 
the particle notion for t >ti [Eq. ((7j)] and the operators 6n(^i), &n(^i) given by 



r 



A+it,)bnit,)+A-it,)blit,) 



(A.2) 



with A^it) defined in (j22I)- Equation (rO|l follows directly from Eq.© with (ET|l . 
(j22)) and ()A.1|) . For motions ending at t = ti with /(ti) ^ Iq the operator bnit) has 
not evolved into the operator An associated with the correct particle notion after the 
dynamics. However, if /(ti) = lo, for example when ti is a multiple of the period of 
boundary vibrations. An = kih)- 



Appendix B. The system for A^'^ and B^f^^ 

Taking the stopping time ti in and (|17|) to be a continuous variable one derives the 
following system of coupled differential equations for Amn and Bmn- 

(B.l) 



Amn 'i^nAmn ^nBmn ~l~ ^ ' \^^nk'^''^^ ^nk^' 

k 

Bmn i^nBmn ~l~ ^nAmn ~l~ ^ ] \^nk^i'nk ^nk'^ 



'mk 



mk 



with 



V (t\ ^^-^^ K^(t\ ^ 



^kit) 

\ ^nit) 



Mnkit) ± 



^nit) 
\ ^kit) 



Mknit) 



(B.2) 



. (B.3) 
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Appendix C. Numerics 

The numerical simulations have been performed by using a Runge-Kutta Prince- 
Dormand method (rkSpd) based on source code provided by the GNU Scientific Library 
(GSL) ^3]. In the table in Fig. K^ll we show the numerical values for the total number 
of particles N{t) created in a cavity subject to sinusoidal oscillations of the form ()28j) 
with frequency u = 3tt [cf. Figs. Hand|21(a)] for two times t = 249.5 and t = 250.0 and 
cut-off parameters /cmax = 20, 30, 40, 50 and 60. The plot in Fig. 1(^11 shows the diagonal 
part of the first of the Bogoliubov relations (fTUj) = 1 — J2m{\^mk\'^ — l-SmfcH = for 
k = 1, ...,10 and 51, ...,60 computed from the solutions of the simulation with A;max = 60. 
The absolute and relative errors for the rkSpd routine in the simulations have been set 
to 10-s. 



X 10 




k 

"'max 


N{t = 249.5) 


N{t -- 


= 250.0) 


20 


0.5799007 


0.51 


^23052 


30 


0.5798943 


0.5t 


^22980 


40 


0.5798951 


0.5t 


^22983 


50 


0.5798956 


0.5i 


^22984 


60 


0.5798959 


0.5i 


^22984 



250 



Figure CI. Left: Plot showing the numerically evaluated diagonal part of the first of 
the BogoHubov relations ifTT))! dk = 1 — X]m(l^™feP~ 1^™^^) for the lowest frequencies 
k — 1,...,10 as well as fc = 51,. ..,60 corresponding to the simulation with cut-off 
^max = 60. Right: Table showing the numerical values of the total particle number at 
times t = 249.5 and t = 250 obtained for w = Stt and fcmax = 20, 30, 40, 50 and 60. 



The plot demonstrates that for those settings rffc = is satisfied by the numerical 
solutions up to ~ 3 X 10~^ at the end of the integration range. Thereby the accuracy 
is better for the lowest modes k = 1, 10 than for the modes k = 51, ...,60. This is 
partly due to the fact that the higher modes are more affected by the truncation of the 
infinite system at /cmax = 60 than the lowest modes. The accuracy for the intermediate 
modes k = 11, 50 lies in between the two "bands" visible in the plot. The remaining 
Bogoliubov relations are satisfied with at least the same accuracy demonstrating that 
the numerical errors are small compared to the values of the particle numbers itself. We 
consider dk as the determining measure for the accuracy of the numerical calculations 
which can be easily enhanced further by increasing the preset accuracy of the integration 
routine. 

The numerical values for N{t = 249.5) and N{t = 250.0) summarized in the 
table in Fig. El] are shown with seven decimal places. Varying the cut-off between 
^max = 30 and 60 both values change only in the last two of the shown seven decimal 
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places and therefore the variation in N{t) when changing /Cmax is smaller than 10~^, 
i.e. smaller than the numerical error in the Bogoliubov relations. This demonstrates 
that the convergence of the numerical values for N{t) is sufficiently good. Furthermore, 
because l{t = 250) = Svre (/q = 1) we can conclude that spurious effects caused by 
discontinuities in the velocity (matching problem) are indeed negligibly small for the 
parameters considered (e <^ 1). 
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